Nonaxisymmetric, multi-region relaxed magnetohydrodynamic equilibrium solutions 
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We describe a magnetohydrodynamic (MHD) constrained energy functional for equilibrium calcu- 
lations that combines the topological constraints of ideal MHD with elements of Taylor relaxation. 
Extremizing states allow for partially chaotic magnetic fields and non-trivial pressure profiles sup- 
ported by a discrete set of ideal interfaces with irrational rotational transforms. Numerical solutions 
are computed using the Stepped Pressure Equilibrium Code, and benchmarks and convergence cal- 
culations are presented. 



I. INTRODUCTION 

The construction of magnetohydrodynamic (MHD) 
equihbria in three-dimensional (3D) configurations is of 
fundamental importance for understanding magnetically 
confined plasmas. To illustrate both the importance and 
subtlety of this problem: it is widely accepted that quies- 
cent plasma confinement depends on constructing equi- 
hbria that are stable to small perturbations, vifhich neces- 
sarily presupposes the existence of an equilibrium; how- 
ever, as pointed out by Grad |l|, "a more primitive reason 
than instability for the lack of confinement is the absence 
of an appropriate equilibrium state" . 

Given that all experimental confinement devices lack 
a continuous symmetry to some extent, either slightly 
so because of discrete coil effects, error fields, or inten- 
tionally applied resonant magnetic perturbations [2|, or 
greatly so because of intrinsic 3D shaping such as in the 
stellarator class 3,^ of confinement devices, the compu- 
tation of 3D equilibrium solutions with arbitrarily chaotic 
fields is of foremost importance. 

The theory and numerical construction of 3D equilib- 
ria is greatly complicated by the fact that toroidal mag- 
netic fields without a continuous symmetry are generally 
a fractal mix of islands, chaotic field lines, and magnetic 
flux surfaces. Any deformation of the plasma bound- 
ary, or error field, that resonates with rational field lines 
will generally (in the absence of ideal surface currents 
at the rational surfaces) result in the formation of mag- 
netic islands. Where these islands overlap, regions of 
connected chaos, so-called stochastic volumes, will form 
0. According to Greene [g], "there is a stochastic re- 
gion in the immediate vicinity of every chain of periodic 
orbits". In contrast, the Kolmogorov-Arnold-Moser the- 
orem indicates that fiux surfaces with "sufficiently irra- 



tional" rotational transform can survive small perturba- 
tion [3, Q. The rotational transform, *, is considered 
sufficiently irrational if it satisfies a Diophantine condi- 
tion, e.g. I i — n/m\ > rm~^ for all integers n and m, 
and where r > and k > 2. 
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FIG. 1: Diophantine pressure profile with r = 0.2 and k = 
2, normalized so that p(0) — 1 and p{l) = 0. The large 
rectangular region is a magnification of the small rectangular 
region. 



The fractal phase-space structure of non-symmetric, 
and therefore generally non-integrable (i.e. such that a 
continuously nested family of fiux surfaces does not ex- 
ist), magnetic fields has important consequences for the 
construction of scalar-pressure ideal equilibria, i.e. so- 
lutions to Vp = j X B. Ideal force balance immediately 
requires that B • Vp = 0, so that the pressure is constant 
along a field line, which in turn implies that the pressure 
must be constant in the stochastic volumes. Pressure 
gradients can be supported on the flux surfaces that sur- 
vive perturbation, but any non-trivial continuous pres- 
sure profile that is consistent with Vp = j x B and the 
fractal structure of a non-integrable field must have a 
gradient that is everywhere discontinuous or zero [9| and 
is something akin to a devil's staircase. For example. 



consider the pressure profile, p{x), defined by p'{x) = 1 
if \x — n/rn\ > rm^^ for all rationals n/m, and p' {x) = 
otherwise. The function p'{x) does not have a well de- 
fined Riemann integral, but an approximation to p{x) is 
shown in Fig. [TJ 

A continuous profile may be constructed, but the gra- 
dient is infinitely discontinuous. This immediately causes 
problems for the construction of scalar-pressure equilib- 
ria, as it is the pressure gradient, Vp, rather than the 
pressure itself, that appears in Vp = j x B. Grad went 
on to conclude that non-trivial equilibria must have "very 
pathological pressure" [1|. The pathological structure of 
scalar-pressure, ideal MHD equilibria with chaotic fields 
causes problems for existing numerical algorithms [yl. 

The most elegant approaches to numerical construc- 
tion of equilibrium solutions employ energy principles, 
and this approach began to be developed as a practical 
method for numerical calculation of 3D ideal-MHD equi- 
libria quite early in the history of computational plasma 
physics [IMl- Today the VMEC code fT^, which has 
been continuously developed since that time, has become 
the most widely used 3D equilibrium code. 

Such codes seek solutions that minimize the plasma 
energy, U = /[p;/(7 — 1) + B^/2^o]dv, but if this func- 
tional is minimized allowing for arbitrary variations with- 
out constraints then the minimizing state is trivial [l5J. 
Instead, in ideal-MHD codes, the plasma energy is min- 
imized under the continuous infinity of constraints im- 
plied by the ideal-MHD equations at each point in the 
plasma, e.g. ^B = V x (^ x B), representing frozen-in 
fiux in each plasma element. However, such variations 
do not allow the topology of magnetic surfaces to change, 
so that reconnection and relaxation phenomena are pre- 
cluded. If the topology assumed in the initial guess for 
the equilibrium is that of a nested family of fiux surfaces 
foliating the plasma domain, then it is constrained to 
remain so during the energy minimization, so that the 
generic state for a 3D equilbrium, one with magnetic is- 
lands and chaos, cannot be obtained in these ideal codes, 
limiting the attainable precision because of the formation 
of unresolved current singularities. 

Thus we need to reduce the number of constraints if we 
are to have a well-posed energy principle for 3D equilib- 
ria, the most extreme reduction being provided by Tay- 
lor relaxation [16]. Allowing for arbitrarily small resis- 
tivity and arbitrary variations, so that the topological 
constraints on the plasma evolution are removed, Tay- 
lor argued that a weakly resistive plasma will relax to 
minimize the plasma energy subject to the constraint of 
conserved helicity, H = J A- 'Bdv. The resulting state is 
a force-free field, and so a globally relaxed Taylor state 
cannot support a pressure gradient — the magnetic field 
topology is not constrained but the pressure is constant 
everywhere, suggesting that a few more constraints are 
required to model a confined plasma satisfactorily. 

Augmenting the helicity constraint with one or more 
smooth moments of A • B leads to a practical approach to 
constructing axisymmetric equilibria [171 Il8| with non- 



trivial pressure profile and low free energy to model a 
toroidally confined plasma with good stability properties. 
However this approach requires good magnetic surfaces 
everywhere in the plasma and is thus limited to inte- 
grable magnetic fields. For the 3D equilibrium problem 
we are forced to consider non-smooth constraints jl9| . 
giving an energy principle that allows for equilibria with 
partially chaotic fields and non-trivial, but discontinuous, 
pressure profile. This model, which we call multi-region, 
relaxed magnetohydrodyamics, or MRXMHD, combines 
ideal topological constraints at a discrete set of selected 
irrational surfaces with Taylor relaxation in between. 
This is a partial relaxation model so that non-zero pres- 
sure gradients are supported at the selected irrational 
surfaces but topological reconnection is possible at the 
intervening rational surfaces. This model is described 
in Sec. nil and in Sec. lIIII we present some illustrations 
of MRXMHD equilibrium solutions, as computed by the 
newly implemented Stepped Pressure Equilibrium Code, 
SPEC. 



II. MULTI-REGION, RELAXED MHD 

The plasma is modeled [13 as a collection of nested an- 
nular regions, Vi for / = 1, .., A'^, which are separated by a 
discrete set of toroidal surfaces, Xi, so that Vi is bounded 
by Ii-i and X;. In each Vj, the magnetic field relaxes 
to minimize the plasma energy, subject to the constraint 
of conserved helicity and mass/entropy. On each Xi , we 
apply the constraints of ideal MHD. Equilibrium states 
are extrema of the constrained energy functional. 



F = Y,{Ui- fiiHi/2 - t^iMi/2) 



(1) 



where the plasma energy, helicity and mass/entropy in 
each annulus are given respectively by 



(2) 
(3) 
(4) 



Ui = 


Jvi \7- 1 2^0/ 


Hi = 


/ A, • B, dv, 


Ml = 


f p]''dv. 



IVi 



The magnetic field is described by a vector potential, 
B; = V X A;, the /i; and ui are Lagrange multipliers, and 
7 is the ratio of specific heats. Hereafter, the permeabil- 
ity of free space factor, /xq, will be omitted. 

In extremizing F we allow for arbitrary variations in 
the pressure, 5pi, and the vector potential, (5A;, in each 
annulus, and the geometry, ^;, of the toroidal surfaces. 
To enforce the constraint that the magnetic field remain 
tangential to the toroidal surfaces, on the I; the varia- 
tions in the field are related to the variations in geometry 
by using a gauge such that 5 A. = ^ x B. We appropri- 
ately call these toroidal surfaces ideal interfaces. 



The first order variation in the free-energy functional 
in each annulus is 
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The equihbrium is thus comprised of a family of nested 
Beltrami fields, where V x B; = /xjB; in each annulus, 
the pressure is constant in each annulus, and the to- 
tal pressure is continuous across the ideal interfaces, 
[[p + B^/2]]=Q. 

Pressure is supported by the ideal interfaces, across 
which a pressure discontinuity is allowed provided there 
is a compensating discontinuity in the tangential field. 
The equilibrium solutions are topologically-constrained, 
partially relaxed, stepped-pressure states. Topological 
constraints on magnetic reconnection have been observed 
in a similar context [201. ^ strong motivation for adopt- 
ing this model is that Bruno & Laurence (2l[ have proved 
that, under certain conditions, such equilibria exist. This 
places the MRXMHD model on a strong mathematical 
foundation. 

In addition to a prescribed boundary, we define 
the equilibrium by prescribing pressure and rotational- 
transform profiles, given as functions of the toroidal flux 
as described below. 

The pressure is constant in each annulus but the rota- 
tional transform, well-defined on flux surfaces and de- 
fined by suitable interpolation across chaotic regions, 
changes across the annuli. Given that the topology of 
the field is arbitrary within the annular regions, and that 
only the interfaces are guaranteed to be flux surfaces, the 
rotational transform is prescribed only on the ideal inter- 
faces. To avoid discontinuous rotational-transform pro- 
files we require the rotational transforms to be the same 
on each side of an interface. So, to specify the pressure 
and rotational-transform profiles, we give the toroidal 
fiux and the pressure, pi, in each annulus, and the ro- 
tational transform, t;, on each interface, for I = l,iV. 
We hold the toroidal flux in each annulus flxed (though 
for reversed-field pinches it would be more appropriate 
to fix the poloidal fluxes). Also, as is typically done in 
equilibrium calculations, we hold the pi and ii constant 
throughout the calculation. However, it should be noted 
that this is not consistent with the variational principle, 
so the constraints Eqs. (jH]) and ^ must be adjusted it- 
eratively during the calculation. 

An initial guess for the geometry of the ideal 
interfaces, I;, is given in cylindrical coordinates, 
{R,<j),Z), via Ri{9,C) ^^-Ri^jCos{mjd — rLjC) and 
Zi{9,Q) — ^- ZijsnY{mj9 ~ TijC,)^ where 9 is at this 
stage an arbitrary poloidal angle and C, = —(p. (We have 
restricted attention to stellarator-symmetric conflgura- 
tions [22].) The toroidal and poloidal fluxes, V't.z and 



^p_i, enclosed by each interface are also assumed given. 
A piecewise-cubic interpolation of the interfaces using 
the radial coordinate s = \p4>t provides a smooth, global 
coordinate system, {s,9,(^), with coordinate Jacobian 
V5 = {Vs-V9x VC)"^ 

The vector potential in each annulus is written 
A — AeW9 + AqWC, and Ag and A^ are discretized us- 
ing a mixed Fourier, finite-element representation, e.g. 
Ag — J2j ^e,j{s) cos{mj9 — rijC), and the radial depen- 
dence is described by A0j[s) —J2k^s,j,k'Pk{s), where 
the ifki-s) are piecewise-quintic basis functions with fi- 
nite support, defined on a radial sub-grid (an example of 
which is shown below). The Agj^k and A^j,fe are con- 
strained to ensure that the flux constraints are satisfled, 
and that ^/gB ■ Vs = dgA^ — d^Ag = at the interfaces, 
but is otherwise general. Setting the derivatives of Eq. ^ 
with respect to the Agj^k and A(^j^k to zero allows each 
B; to be efficiently determined as the solution to a sparse 
system of linear equations. Each B/ depends only on the 
geometry of bounding interfaces, J;_i and Ii, the en- 
closed fluxes, and the Lagrange multiplier, /zj, which is 
related to the parallel current. This must be adjusted 
to ensure that the rotational transform on the ideal in- 
terfaces matches the prescribed value. The computation 
of the Beltrami fields in multiple regions is trivially dis- 
tributed across multiple processors. 

The innermost volume contains the coordinate origin, 
where the coordinate Jacobian goes to zero. At the ori- 
gin, we enforce the condition that the geometry of the 
interfaces is regular, and the geometry of the innermost 
interface is obtained by extrapolation. 

The interface geometry must be adjusted in order to 
satisfy force balance, [[p + B'^/2]] = 0. The first-order 
variation in the energy functional, F, depends only the 
normal component of the geometrical variation, ^ • n. 
In order to obtain a unique Fourier representation of 
the interface geometry, we follow the approach used in 
VMEC [231 and construct an angle that minimizes a 
measure of the spectral width, Y^AmFj + n'j){R'j + Zj), 
where p and q are arbitrary integers controlhng the 
degree of spectral condensation (in the following section 
we choose p = 4 and g = 4), and so obtain an opti- 
mally accurate representation of the interface geometry 
with a finite set of Fourier harmonics. Allowing for 
tangential variations of the form, SR — dgRSu and 
6Z — dgZ Su, the condition that 9 minimizes the spec- 
tral width of each interface is // = dgRi X + dgZiY = 0, 
where X = ^ (m^ -|- nj)i?;j cos(mj0 — n^C) and 
Y = E/"^j + n'j)Zij sm{mj9 - UjO- 

The interface geometries are described by the Rij and 
the Zij, which we collect together as a vector, x. We con- 
struct a vector of constraints, JF, as a collection of the 
Fourier harmonics of the force imbalance, [[p-|-i3^/2]]i 



\i,j^ 



and the spectral constraints, Iij, at each interface. The 
task of constructing equilibrium solutions is thus reduced 
to the standard mathematical problem of finding a zero of 
a multi-dimensional function, ^(x) = 0, which is solved 
using a mixed Newton, convex-gradient method provided 



by the NAG library. Further details of the algorithm, in- 
cluding convergence studies and benchmark calculations, 
will be presented elsewhere. 

An analysis of the force-balance condition, 
[[p + B"^ /2]]=Q, shows that, generally, in order for 
an interface to support pressure, it must have irrational 
rotational transforms on each side [24]. As previously 
mentioned, we take these to be equal. "Noble" irra- 
tionals play an important role in chaos as, typically, 
flux surfaces with noble rotational transform are lo- 
cally the most robust [y, [2^. We thus constrain the 
ideal interfaces that support pressure to have noble 
rotational transform, given as a Fibonacci-ratio limit: 

Pn+l/qn+l = [Pn-l + Pn) / {Qn-l + Qn) aS n -> OO, 

beginning from any pair of rationals, pi/qi and P2I12 
which satisfy \p\qi—P2(li\ = 1- The limiting ratio is 
* = (Pi + lP2)/{qi + 792), where 7 = (1 -I- •\/5)/2 is the 
golden mean. As mentioned earlier, adjustments must 
be made iteratively during the equilibration calculation. 
Specifically, we adjust the helicity multiplier, /i/, which 
is related to the parallel current, and the poloidal flux in 
each annulus to satisfy the interface rotational-transform 
constraints. 



III. NON-AXISYMMETRIC EXAMPLES 




FIG. 2: Smooth pressure profile, p{ip) ~ 1 — 2i/) + 'i/'^, supplied 
to VMEC, and the stepped pressure profile used in SPEC. 




In the following we present numerical examples of 
stepped-pressure equilibria with 3D boundaries. The 
following examples are constructed by applying a non- 
axisymmetric deformation to an otherwise axisymmet- 
ric configuration. The axisymmetric configuration is de- 
fined by a fixed boundary with a major radius of 1 m 
and a circular cross section with minor radius 30 cm. For 
simplicity of illustration, and computational expediency, 
we restrict attention to equilibria with only four annular 
regions, bounded by four interfaces. The pressure pro- 
file, shown in Fig. [51 is a piecewise-flat approximation 
to p(^) = po(l — 2?/; -I- 1/)^), where po is a scaling factor. 
The interface rotational-transform profile is a discrete ap- 
proximation to i{il)) =0.8839642543 -0.7799929021V', as 
described in Table. U and shown as the small squares in 
Fig. 11 

TABLE I: Ideal-interface rotational transform 



0.0 



0.2 



0.4 



0.6 



0.8 



1.0 





f 


V^t 
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(6 + 77)/(7-f78) 


= 0.8687325 . . . 


0.0195280 
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(2 + 73)/(3 + 74) 


= 0.7236068 . . . 


0.2055884 


3) 


(l + 7l)/(2 + 73) 


= 0.3819660 . . . 


0.6435933 


4) 


(l+7l)/(9 + 7lO) 


= 0.1039714 . . . 


1.0000000 



Below we examine two non-axisymmetric configura- 
tions. The first is a zero-beta configuration (i.e. po — 0) 
with an TO = 0, n = 1 deformation of the boundary. This 
results in an equilibrium with a helical magnetic axis. In 
the second "high-pressure" equilibrium, po is increased 



FIG. 3: rotational-transform profile (black line) for the ax- 
isymmetric "high-pressure" configuration shown in Fig. |4l 
The small squares indicate where the rotational transform 
is constrained. The smooth profile (grey line) is supplied to 
VMEC. 



in order to induce a 3 cm Shafranov shift of the magnetic 
axis of the axisymmetric configuration. 

At this point we stress that, despite the fact that only 
4 interfaces were used and the piecewise-flat approxima- 
tion to the smooth pressure profile seems crude, this is 
sufficient to describe the effect of the pressure on the 
global geometry of the equilibrium. To illustrate this, 
we compare the axisymmetric high-pressure equilibrium 
computed with SPEC (using the stepped pressure profile) 
to the corresponding VMEC equilibrium (computed with 
the smooth pressure profile) . For the VMEC calculation, 
the smooth rotational-transform profile shown in Fig. [3] 
was used. Shown in Fig. |4] are the cross-sections of the 
interfaces as computed by SPEC and the corresponding 
irrational flux surfaces computed by VMEC. That the ge- 
ometry of the interfaces agrees so well is partly due to the 
fact that the Shafranov shift is rather insensitive to the 
pressure profile itself, but depends on the integral of the 
pressure, i.e. the plasma beta. Also, the location of the 
magnetic axis and the geometry of the innermost inter- 
face computed by SPEC agree well with that computed 



by VMEC. 

In Fig. 13] we also show the radial sub-grid resolution 
used in each annulus. In total, there are 78 global radial 
degrees of freedom used in the piecewise-quintic represen- 
tation of each Fourier harmonic of the magnetic vector 
potential. 

Between the interfaces in the SPEC equilibrium, the 
rotational transform is not given a priori: the rotational 
transform within each annulus is to be determined as 
part of the equilibrium calculation. An approximation 
to the global rotational transform may be constructed 
a posteori by field-line following, and this is shown in 
Fig.El 

To this high-pressure axisymmetric equilibrium we add 
a boundary perturbation that resonates with low-order 
rational surfaces inside the plasma, which induces a large 
fraction of the equilibrium magnetic field to become 
chaotic. While the size of the magnetic islands created 
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FIG. 4: The ideal interfaces of the high-pressure equilibrium 
computed by SPEC (thin lines, both upper and lower half) 
and the corresponding irrational surfaces of the VMEC equi- 
librium (thick lines, upper half only). In the lower half, the 
radial sub-grid and the computational angle coordinate is also 
shown (grey lines). 

by 3D shaping is determined by the magnitude of the res- 
onant component of the deformation, and the shear, their 
location is determined by the rotational-transform pro- 
file. The rotational transform is specified discretely and 
is only constrained at the interfaces themselves; however, 
with the interface rotational transforms given in Table. U 
the location of the t = 1/2 resonance is guaranteed to lie 
in the third annulus (that lies between the second and 
third interfaces). Other n ~ 1 resonances, such as the 
i = 1/3, i = 1/4 and * — 1/5 resonances, are guaran- 
teed to lie in the fourth (outermost) annulus. The m = 0, 
n = 1 resonance is not present. 

In all of the calculations referred to below, force bal- 
ance at each interface is satisfied to within IjFI < 10^^^, 



or less. 



A. Helical magnetic axis 

For the zero-beta equilibrium, with po — 0, the fixed 
outer boundary is described by 

R = 1.0 + 5cos(C) + 0.3cos(6i), 
Z = (5sin(C)-l-0.3sin(6l), 

where the helical deformation is (5 = 0.035. The Fourier 
representation of both the ideal interfaces and the vector 
potential in each annulus includes the modes < m < M 
and -N <n< N, for {M,N) = (6,3). Poincare plots 
of the equilibrium are shown on three cross sections in 

Fig.m 

Nonaxisymmetric configurations are not guaranteed to 
be integrable, but neither are they guaranteed to be glob- 
ally chaotic. If one were to look closely, small islands may 
be observed at all the rational surfaces inside the plasma; 
but in this case, because the nonaxisymmetric deforma- 
tion of the boundary does not directly resonate with any 
rational surfaces, no large island chains will form. 

To illustrate an equilibrium that does have a signif- 
icant volume of chaotic fields, in the next example we 
include large boundary deformations that resonate with 
low order rational surfaces. 



B. Strongly chaotic equilibrium 

To drive islands in the high-pressure configuration we 
include a deformation of the minor radius that resonates 
directly with the t = 1/2 and i = 1/3 rational surfaces. 
The outer boundary is given by 



R = 
Z = 



1.0 



[0.3 + S cos{20 - 
[0.3 + (5cos(26» 



-C) 



Scos{36 - 
' dcos{39 



C)]cos{e), 



where the magnitude of the deformation is given 
S = 0.003. The Fourier representation includes the modes 
< m < M and -iV < n < A^, for (M, A^) = (9, 4). The 
field in the outermost annulus is now strongly chaotic. 
This is because this annulus contains several low-order 
resonances, e.g. t = 1/3, which is directly driven by 
the applied boundary deformation, and also the t = 1/4, 
t = 1/5, and i — 1/6, and islands will form at these lo- 
cations islands because of toroidal and poloidal coupling. 
These islands are quite close together and the magnitude 
of the deformation is sufficient to ensure that these is- 
lands overlap. 

To confirm that this solution is converged with respect 
to Fourier resolution, we recompute this equilibrium us- 
ing the Fourier resolution {M,N) = (6,1), (7,2), and 
(8,3), and compare the geometry of the interior inter- 
faces. Simply comparing the Fourier harmonics of the 
interfaces at different resolutions may give misleading re- 
sults because, as the Fourier resolution is increased, the 
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FIG. 5: Poincare plots (gray dots) of the helical axis equilibrium on (A) ( — 0; (B) ( — n/2; and (C) (^ — n. The ideal interfaces 
are shown with black lines. 





FIG. 6; Poincare plots (gray dots) of the resonantly deformed equilibrium on the (A) ( — 0, (B) ( — n/2, and (C) ( — n 
toroidal cross-sections. The ideal interfaces are shown with black lines. 



spectral condensation algorithm has more opportunity to 
exploit the tangential freedom, and thus to give a slightly 
different angle parametrization of the same geometrical 
interface. To eliminate any uncertainty arising from this, 
we introduce the following angle-independent measure of 
the geometrical difference between a given pair of inter- 
faces. 



dM,N{d,a) = y[x{0) - XM,N{a)f + [y{0) - yM,N{a)] 



of the error according to 
r2v 



^AI,N 



JQ 



(6) 



where Dm,n{0) = miuQ dM,N{0, oi). Taking the reference 
configuration to be the solution computed with the high- 
est Fourier resolution, i.e. with {M,N) = (9,4), so that 
^9.4 = by definition, the convergence error for the in- 
terior interfaces is shown in Fig. [T] 



where x{9) = R{9Xo) and y{9) = Z{9,Cq) is the inter- 
face cross-section curve of a reference solution (specified 
below), and XM,Nia) = -Rj\/,Ar(a, Co) and yM,N{a) = 
Zm,n{o!^,Co) is the cross-section curve of the solution 
computed with Fourier resolution {M,N)^ on the plane 
C = Co ■ We then compute the angle-independent measure 



IV. COMMENTS 

In this article we have constructed stepped-pressure 
equilibria with only 4 nested annular regions, as this 
seems to be sufficient in order to capture the global de- 
formation induced by non-trivial pressure, as confirmed 
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FIG. 7: Angle-independent measure of the error in the interior 
interface cross-section geometry, with reference configuration 
(M, iV) = (9, 4) and Co = 0. 



by Fig. m It is possible to compute equilibria with 
arbitrarily many interfaces and annular regions; how- 
ever, as the number of ideal interfaces increases, and 
the separation between the ideal interfaces decreases, the 
present modified-Newton algorithm becomes more fragile 
because of the emergence of small eigenvalues in the ma- 
trix, V^. In future work we hope to explore alternative 
numerical algorithms for constructing solutions. Given 
that stepped-pressure equilibria are defined as extrema 
of a constrained energy functional, it should be possi- 
ble to implement a rapid, preconditioned descent-style 
algorithm |2q . Also in future work, we intend to com- 
pute the Beltrami field in the innermost volume (which 
contains the toroidal coordinate singularity), so that the 
geometry of the innermost interface can be determined 
directly from force balance, rather than appealing to reg- 
ularity conditions at the origin and employing extrapo- 
lation methods. 

Specifying the profiles discretely may seem arbitrary, 
but it is a practical means of maintaining some a priori 
control over the pressure and rotational transform while 
making minimal assumptions regarding the topology of 
the field. In fact, the only topological assumption made 
is that, where there are pressure gradients there must 
be irrational flux surfaces. Thus the prescription of the 
profiles may be made arbitrarily fine-grained as long as 
magnetic surfaces exist. Where they do not exist, force- 
free fields are the only choice consistent with no mass 
flow, Beltrami fields being the most practical. 

To conclude, we make some comments regarding the 
existence of MHD equilibrium solutions, and the "patho- 
logical" nature of the pressure profile. Earlier we com- 
mented how solutions to Vp = j x B with non-integrable 
magnetic fields with a fractal phase space must have pres- 
sure profiles with infinitely discontinuous gradient. In 
the MRXMHD model, the pressure profile is piecewise 
flat, and possibly discontinuous at the ideal interfaces. 



Such a profile may also be described as pathological; 
however, the MRXMHD model is based on an integral 
principle, and a discontinuous pressure profile remains 
an integrable function, provide the number of discontinu- 
ities is finite, as it certainly is in the above calculations. 
The MRXMHD equilibrium is well defined mathemati- 
cally, and at no point in the numerical construction of 
the stepped-pressure equilibrium is the pressure gradient 
required. 

In the analysis of the force-balance condition, 
[[p + -B^/2]] = 0, arising in the Euler-Lagrange equation, 
Eq. ([5]) , it was shown [2J| that generally pressure can only 
be supported if the interfaces have irrational rotational 
transform. This in turn places constraints on the pres- 
sure and rotational-transform profiles that are used to 
define the equilibrium: if pressure is placed on the ratio- 
nal interfaces, then no equilibrium solution will exist. 

An analogous condition holds for ideal, scalar pres- 
sure equilibria with nested fiux surfaces, i.e. integrable 
magnetic fields. States that minimize the plasma en- 
ergy, U, allowing for ideal variations, must satisfy the 
Euler-Lagrange equation, Wp = j x B, which is the anal- 
ogous statement of force-balance for ideal MHD equilib- 
ria with nested flux surfaces, i.e. ideal equilibria that 
are globally topologically constrained, rather than dis- 
cretely topologically constrained as in MRXMHD. An 
analysis of this equation shows that there is a singular- 
ity in the resonant harmonic of the parallel current at 
each rational surface j27ij , which are dense in any system 
with shear. In addition to the ideal ^-function surface 
currents required to shield resonant flelds, that would 
otherwise result in the formation of magnetic islands, 
there generally exist pressure-driven 1/x style singulari- 
ties, where x = (t — n/m), which are required to satisfy 
quasi-neutrality. Writing the current as j = crB + jj_, 
and insisting that V • j = 0, the parallel current, a, must 
satisfy B • Vct = —V • j_l, where the perpendicular cur- 
rent is given by force balance. The magnetic differential 
equation is singular, and solvability conditions on the 
perpendicular current must be satisfled if a single val- 
ued (7 is to exist |28j : an arbitrary j^ = B x Vp/B^ 
is not consistent with quasi-neutrality! The singularity 
in the B • V operator is exposed by employing straight- 
fleld-line coordinates, which can be constructed globally 
if, and only if, the magnetic held is integrable, so that 
Y^B -V = 8^+ ide ■ The solvability condition that must 
be satisfied for quasineutrality is that, in arbitrary ge- 
ometry, the pressure gradient must go to zero at each 
rational surface at least as fast as ( a — n/m) , and at least 
as fast as ( i — n/m)'^ if the pressure is to remain mono- 
tonic. Given that the rational surfaces are dense (i.e. 
arbitrarly close to any point in space) in plasma equilib- 
ria with shear, this results in a pressure profile that may 
also be described as pathological. 

The MRXMHD approach seeks integrable solutions, 
rather than differentiable solutions. The philosophy of 
seeking weak solutions was endorsed by Garabedian, who 
claimed that "differentiable solutions of the equilibrium 



equations do not exist in general when the geometry is 
three-dimensional, so that weak solutions are required to 



model the physics adequately" [29| . 
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